Lecture 7#
Interpolation#
We have up until now focussed solely on how to solve equations using finite difference stencils and mesh functions evaluated on a regular (uniform) mesh. For example, in 1D a mesh function is
where \(x_i = i \Delta x\). In 2D a mesh function is
where \(y_j = j \Delta y\). We have also looked at a problem that is 1D in space, plus time. In this case the mesh function is
where \(t_n = n \Delta t\). Later in this lecture we will look at the wave equation in 2D space plus time. In this case the mesh function will be
Up until now we have learned how to compute these mesh functions, but we have not really talked about how to use them. A natural question is what if you need to know the solution at a specific point and time \(u(t, x, y)\)? And this point is not a mesh point? How do you compute this number from the mesh function?
The answer is interpolation.
One-dimensional interpolation#
Interpolation is most easily understood using a one-dimensional example. Assume that the mesh is \(\Omega \in [0, 1]\) and
Let us create a mesh using \(N=5\) and the mesh function corresponding to (7). Let us also plot the mesh functions as blue dots.
Show code cell source
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
N = 5
xj = np.linspace(0, 1, N+1)
u = xj*(1-xj)
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 3))
ax0.plot(xj, u, 'bo')
ax1.plot(xj, u, 'bo-');
In the figure above to the left the mesh function \(U=(u_i)_{i=0}^N\) is represented with the blue dots. In the figure to the right linear profiles have been added by matplotlib in between the points because we ask for it using bo-. The linear profiles represent linear interpolation, which is a first order approach to computing the mesh function in between the mesh points. The next figure compares the linear interpolation to the exact solution.
fig = plt.figure(figsize=(5, 3))
plt.plot(xj, u, 'bo-')
xl = np.linspace(0, 1, 1000)
plt.plot(xl, xl*(1-xl), 'r-');
Note that even though the mesh function here is exact at all the mesh points, the linear interpolation between points is far from being accurate.
We now want to compute \(u(0.75)\) using linear interpolation. We get this from drawing a straight line from \(u(x_3=0.6) = 0.24\) to \(u(x_4=0.8)=0.16\), and this line can be represented as
By inserting for \(u_3, u_4, x_3\) and \(x_4\) we get
We can create this profile using Sympy
x = sp.Symbol('x')
uo = u[3] + (u[4]-u[3])/(xj[4]-xj[3])*(x-xj[3])
uo
And we can plot the linearly interpolated point \(u(0.75)\) as a green dot
fig = plt.figure(figsize=(5, 3))
plt.plot(xj, u, 'bo-')
xl = np.linspace(0, 1, 1000)
plt.plot(xl, xl*(1-xl), 'r-')
plt.plot(0.75, uo.subs(x, 0.75), 'go');
The interpolation curve that we have created using Sympy can be used for any point, but it makes most sense to use it between \(x_3\) and \(x_4\). If we use it outside this domain, the accuracy is much worse, and we also refer to this as extrapolation.
fig = plt.figure(figsize=(5, 3))
plt.plot(xj, u, 'bo-')
xl = np.linspace(0, 1, 1000)
plt.plot(xl, xl*(1-xl), 'r-')
plt.plot(xl[500:], sp.lambdify(x, uo)(xl[500:]), 'g-');
Lagrange interpolation polynomials#
Linear interpolation is easy to understand and easy to compute with. However, it is not very accurate. And why would you want to waste your hard earned second-order accurate mesh function (or exact, as above) using only first order accurate interpolation? You shouldn’t and there is no need to. Linear interpolation makes use of only 2 interpolation points. Higher order interpolation simply makes use of more points. The most straight-forward way to do this is through the use of Lagrange interpolation polynomials.
A Lagrange interpolation polynomial makes use of \(k+1\) nodes \((x^0, x^1, \ldots, x^k)\). These are not usually identical to the entire computational mesh, but rather the closest points to the point of interest. Like for the linear interpolation above, the points would simply be \((x^0, x^1) = (x_3, x_4)\) for linear interpolation, or \((x^0, x^1, x^2) = (x_2, x_3, x_4)\) for second order.
Note
To separate the Lagrange points from the mesh points we will simply use superscript on the Lagrange points. So here the superscript does not mean the power!
With Lagrange interpolation we make use of the Lagrange basis functions (or cardinal functions)
which can be more compactly written as
The Lagrange basis functions can be easily implemented using Sympy:
def Lagrangebasis(xj, x=x):
"""Construct Lagrange basis function for points in xj
Parameters
----------
xj : array
Interpolation points
x : Sympy Symbol
Returns
-------
Lagrange basis functions
"""
from sympy import Mul
n = len(xj)
ell = []
numert = Mul(*[x - xj[i] for i in range(n)])
for i in range(n):
numer = numert/(x - xj[i])
denom = Mul(*[(xj[i] - xj[j]) for j in range(n) if i != j])
ell.append(numer/denom)
return ell
We can plot some of the basis functions, for example for second order interpolation using \((x^0, x^1, x^2) = (x_2, x_3, x_4)\)
ell = Lagrangebasis(xj[2:5], x=x)
for l in ell:
print(l)
12.5*(x - 0.8)*(x - 0.6)
-25.0*(x - 0.8)*(x - 0.4)
12.5*(x - 0.6)*(x - 0.4)
fig = plt.figure(figsize=(5, 3))
plt.plot(xj, u, 'bo-')
xl = np.linspace(xj[2], xj[4], 100)
plt.plot(xl, sp.lambdify(x, ell[0])(xl), 'r:')
plt.plot(xl, sp.lambdify(x, ell[1])(xl), 'g:')
plt.plot(xl, sp.lambdify(x, ell[2])(xl), 'b:')
plt.legend(['u(x_i)', '$\ell_0(x)$', '$\ell_1(x)$', '$\ell_2(x)$'], loc='upper left');
Note
Note that all Lagrange basis functions are such that on the chosen mesh points \((x^{i})_{i=0}^{k}\) we have
Note
The Lagrange basis functions do not depend on the mesh function values, only on the mesh points. Also note that the mesh points do not need to be uniform, any mesh will do as long as all points are different, i.e., \(x^k \ne x^j\) for \(k \ne j\)
The Lagrange interpolating polynomial is defined as
where \(u^j\) are the mesh function values at the chosen mesh points. Note that for the linear interpolation example, we would use the two neighbouring values \((u^0, u^1) = (u_3, u_4)\).
Because of (8) we automatically get that
since
and \(\ell_j(x^i) = 1\) for \(i=j\), whereas \(\ell_j(x^i) = 0\) for all \(i \ne j\).
The Python function Lagrangefunction returns the complete Lagrange function:
def Lagrangefunction(u, basis):
"""Return Lagrange polynomial
Parameters
----------
u : array
Mesh function values
basis : tuple of Lagrange basis functions
Output from Lagrangebasis
"""
f = 0
for j, uj in enumerate(u):
f += basis[j]*uj
return f
Let us try this new function and see if we can reproduce the linear interpolation function that we have already computed above:
ell = Lagrangebasis(xj[3:5])
L = Lagrangefunction(u[3:5], ell)
print(L)
0.48 - 0.4*x
Works!
Term by term this is
The Lagrange polynomial is then
Inserting for \(u^0=0.24, u^1=0.16, x^0=0.6\) and \(x^1=0.8\) we get
which can be rearranged into
which is the same as the linear interpolation we have already computed before.
Note
Higher order interpolation simply makes use of more neighbouring points.
ell = Lagrangebasis(xj[2:6])
L = Lagrangefunction(u[2:6], ell)
L
We can evaluate the Lagrange polynomial for an array of points as well:
sp.lambdify(x, L)(np.array([0.5, 0.6, 0.7]))
array([0.25, 0.24, 0.21])
Two-dimensional interpolation#
We now consider the 2D problem where we have the mesh function
and want to evaluate this mesh function for some point inside the domain that is not a mesh point. We now need to interpolate in two dimensions. To get started we need a 2D mesh and an appropriate mesh function to interpolate. We will use the 2D function
and a domain \(\Omega = [0, 1]^2\).
Note
It is common to write \([a, b]^2\) for the Cartesian product of a set with itself. So \([0, 1]^2 = [0, 1] \times [0, 1]\)
def mesh2D(Nx, Ny, Lx, Ly, sparse=False):
x = np.linspace(0, Lx, Nx+1)
y = np.linspace(0, Ly, Ny+1)
return np.meshgrid(x, y, indexing='ij', sparse=sparse)
N = 10
xij, yij = mesh2D(N, N, 1, 1, False)
u2 = xij*(1-xij)*yij*(1-yij)
plt.figure(figsize=(3, 2))
plt.contourf(xij, yij, u2);
We want to find \(u(0.55, 0.65)\) from the mesh function \(U=(u_{ij})_{i,j=0}^{N, N}\).
There are several procedures for performing 2D interpolation. We will choose the Lagrange interpolation first and use \(k+1\) points in both \(x\) and \(y\) directions. In this case the Lagrange interpolating polynomial will be
where \(u^{m,n}\) are the mesh function values for the points we choose for the interpolation, that should be the points closest to the interpolated point. For our case with \(N=10\) this would be the four points \((u_{ij})_{i=5, j=6}^{6, 7}\). Using only these four points we will achieve linear interpolation in both \(x\) and \(y\)-direction.
The four interpolation points are shown as blue dots in the figure below, and the point we are trying to evaluate is shown as a red dot
Show code cell source
plt.figure(figsize=(6, 4))
plt.contourf(xij, yij, u2)
plt.scatter(xij[5:7, 6:8], yij[5:7, 6:8])
plt.plot(np.array([0.55]), np.array([0.65]), 'ro')
plt.xticks(xij[:, 0])
plt.yticks(xij[:, 0])
plt.grid()
The Lagrange basis functions \(\ell_m(x)\) and \(\ell_n(y)\) are computed exactly as before, only now with \(\ell_m(x)\) using points along the \(x\)-direction (for linear \((x^0, x^1)=(x_5, x_6))\) whereas \(\ell_n(y)\) are using points along the \(y\)-direction (for linear \((y^0, y^1)=(y_6, y_7)\)). Naturally, we can choose more points if we want more than linear interpolation.
We can reuse the previously defined function Lagrangebasis, since both \(\ell_m(x)\) and \(\ell_n(y)\) are one-dimensional functions. This was the reason why we gave Lagrangebasis a keyword argument x=x, because now we can compute \(\ell_n(y)\) using a Symbol y as shown below:
y = sp.Symbol('y')
lx = Lagrangebasis(xij[5:7, 0], x=x)
ly = Lagrangebasis(yij[0, 6:8], x=y)
We need a new function for the 2D Lagrange polynomial though, since this is a double loop
def Lagrangefunction2D(u, basisx, basisy):
N, M = u.shape
f = 0
for i in range(N):
for j in range(M):
f += basisx[i]*basisy[j]*u[i, j]
return f
Lagrangefunction2D should return a Lagrange function like (9).
We can now test the function by first creating it
f = Lagrangefunction2D(u2[5:7, 6:8], lx, ly)
f
sp.simplify(f)
and then checking that it works
ue = x*(1-x)*y*(1-y)
print(f.subs({x: 0.55, y: 0.65}), ue.subs({x: 0.55, y: 0.65}))
0.0551250000000000 0.0563062500000000
The results are not identical because this is merely a first order interpolation. With one additional point using \((x_5, x_6, x_7)\) and \((y_5, y_6, y_7)\) we get
lx = Lagrangebasis(xij[5:8, 0], x=x)
ly = Lagrangebasis(yij[0, 5:8], x=y)
f = Lagrangefunction2D(u2[5:8, 5:8], lx, ly)
print(f.subs({x: 0.55, y: 0.65}), ue.subs({x: 0.55, y: 0.65}))
0.0563062500000000 0.0563062500000000
These two results are now identical since the polynomials are second order.
Other tools for interpolation#
There are also other tools in Numpy and Scipy that can help you do 2D interpolation with minimal effort. Most of them are based on splines, which is not a topic in this class, but still useful to know about. We can do interpolation with simple, one-line calls, as follows
from scipy.interpolate import interpn
print(interpn((xij[5:8, 0], yij[0, 5:8]), u2[5:8, 5:8], np.array([0.55, 0.65])))
[0.055125]
Here we get the same result as linear above even though we use 3 points in each direction. This is because interpn is using linear interpolation by default. Modify this to cubic (which requires 4 points in each direction) to obtain better accuracy
print(interpn((xij[5:9, 0], yij[0, 5:9]), u2[5:9, 5:9], np.array([0.55, 0.65]), method='cubic'))
[0.05630625]
Errors and integration#
For the numerical solution \(u\) and the exact solution \(u^{e}\) the \(L^2\) error norm can in general, for any domain \(\Omega\) and number of dimensions, be defined as
In 2D \(\Omega = [0, L_x] \times [0, L_y]\) and the integral becomes
Hence we need to be able to compute such double integrals. One approach is the midpoint rule computed for any integrand \(u(x, y)\) as follows
This integral requires the value of the function \(u\) in the middle of all the computational cells, where the cells are defined as all the squares in the figure below. In order to compute the value of \(u\) in the center of all cells, we need to interpolate. Linear interpolation will in this case lead to an average of \(u\) in the four surrounding nodes, marked with blue dots in the figure for \(u(2.5 \Delta x, 1.5 \Delta y)\).
Show code cell source
fig = plt.figure(figsize=(4, 4))
N = 5
for i in range(N):
plt.plot((0, N-1), (i, i), 'k')
for j in range(N):
plt.plot((j, j), (0, N-1), 'k')
ax = fig.gca()
ax.axis('off')
ax.axis('equal')
plt.plot(2.5, 1.5, 'ko')
plt.plot(2, 1, 'bo')
plt.plot(2, 2, 'bo')
plt.plot(3, 1, 'bo')
plt.plot(3, 2, 'bo')
plt.text(2.5, 1.6, r"$u(\frac{5\Delta x}{2}, \frac{3\Delta y}{2})$", horizontalalignment='center');
Lets create a mesh function from an analytical solution
N = 20
xij, yij = mesh2D(N, N, 1, 1, False)
u2 = np.cos(xij)*(1-xij)*np.sin(yij)*(1-yij)
ue = sp.cos(x)*(1-x)*sp.sin(y)*(1-y)
We can compute values of \(u\) in all the cell centers easily using vectorization. For the field u2 defined above with shape \((N_x+1)\times (N_y+1)\) we get
um = (u2[:-1, :-1]+u2[1:, :-1]+u2[:-1, 1:]+u2[1:, 1:])/4
Note that um is a matrix of shape \(N_x \times N_y\).
The integral required for the \(L^2\)-norm can now be computed easily and we get the following function for the \(L^2\) error
def L2_error(u, dx, dy):
um = (u[:-1, :-1]+u[1:, :-1]+u[:-1, 1:]+u[1:, 1:])/4
return np.sqrt(np.sum(um**2*dx*dy))
ua = sp.lambdify((x, y), ue)(xij, yij)
L2_error(u2-ua, 1/N, 1/N)
4.961003963114361e-18
The result is close to machine precision, which it should be because the Lagrange interpolator produces results identical to the mesh function at all mesh points. Hence, to test the integral, we can send in just the mesh function instead:
L2_error(u2, 1/N, 1/N)
0.09552823074589575
The exact answer is close
float(sp.sqrt(sp.integrate(sp.integrate(ue**2, (y, 0, 1)), (x, 0, 1))))
0.09586332023451888
Just like for the 1D case there is a simplification using the small \(\ell^2\) norm instead
The result will be very similar to the \(L^2\)-norm for most mesh functions. Note how the sums are from 0 to \(N_x\) and \(N_y\) above, since these are sums over nodes and not cell centers. The \(l^2\) error is very easy to implement:
def l2_error(u, dx, dy):
return np.sqrt(dx*dy*np.sum(u**2))
l2_error(u2, 1/N, 1/N)
0.09980275105614006
which is less accurate than the midpoint rule, but still useful as an error measure.
You can also use nested trapz or Simpson’s rule for the integration:
def L2_error_trapz(u, dx, dy):
return np.sqrt(np.trapz(np.trapz(u**2, dx=dy, axis=1), dx=dx))
L2_error_trapz(u2, 1/N, 1/N)
0.09592899344305951
def L2_error_simps(u, dx, dy):
from scipy.integrate import simpson as simp
return np.sqrt(simp(simp(u**2, dx=dy, axis=1), dx=dx))
L2_error_simps(u2, 1/N, 1/N)
0.09586418448463656
The wave equation in 2D plus time#
We will now consider the wave equation
in a two-dimensional spatial domain \(\Omega = [0, L_x] \times [0, L_y]\) and in time \(t \in [0, T]\). The problem will be solved with Dirichlet boundary conditions \(u=0\) on the entire boundary, and it will also require an initial condition \(u(0, x, y) = I(x, y)\) and an initial derivative, here set to zero \(\frac{\partial u}{\partial t}(0, x, y) = 0\).
Perhaps surprisingly, this equation is considerably easier to implement than the Poisson equation, because of the hyperbolic nature of the problem. Hyperbolic problems are often solved with recursive methods that simply step the solution forward in time, and do as such not require setting up and solving matrix equations. However, we can still use much of the theory described in lecture 6, about the 2D mesh and 2D discretised equations.
Note
Time-dependent equations like the wave-equation (hyperbolic) are usually easy to implement, but the numerical scheme needs to consider stability as the solution is moved forward in time. Steady-state equations like the Poisson equation (elliptic) are usually more difficult to set up, but stability is then usually not an issue because we are not stepping the solution anywhere (unless we try to solve the problem iteratively).
Mesh function#
The computational mesh is the same as in lecture 6, with the domain
where \(\boldsymbol{x} = (x_i)_{i=0}^{N_x}\), \(\boldsymbol{y} = (y_j)_{j=0}^{N_y}\),
and \(\Delta x = L_x/N_x\) and \(\Delta y = L_y/N_y\).
Time will also use uniform intervals, and \(t_n = n \Delta t\) for \(n=0, 1, \ldots, N_t\). A function \(u(t, x, y)\) evaluated on this mesh is now found as
However, we will only store the solution at three consecutive time steps, just like for the 1D wave problem. We will as such make use of the three mesh functions
that all are matrices of shape \((N_x+1) \times (N_y+1)\).
We will solve the wave equation using central difference in time as well as in space. A discretization of all internal points is thus
On matrix form this equation reads
The initial condition \(\frac{\partial u}{\partial t}(0, x, y) = 0\) can be implemented like we did in lecture 5 for the wave equation with one spatial dimension. We use
in the PDE for \(n=0\)
such that
A solution algorithm is thus to
Specify \(U^0\) and \(U^1\) from initial conditions
for n in (1, 2, …, \(N_t-1\)) compute
\(U^{n+1} = 2U^n - U^{n-1} + (c\Delta t)^2 \left( D^{(2)}_x U^n + U^n (D^{(2)}_y)^T \right)\)
Apply boundary conditions to \(U^{n+1}\)
Swap \(U^{n-1} \leftarrow U^n\) and \(U^n \leftarrow U^{n+1}\)
Other than that we can simply reuse a lot of code from lecture 6. A function for the second derivative matrix is:
import numpy as np
import scipy.sparse as sparse
def D2(N):
D = sparse.diags([1, -2, 1], [-1, 0, 1], (N+1, N+1), 'lil')
D[0, :4] = 2, -5, 4, -1
D[-1, -4:] = -1, 4, -5, 2
return D
A function mesh2D for creating the 2D mesh has already been described in this notebook. We are thus ready to create a solver that steps the solution forward in time. For simplicity we will use \(N_x=N_y=N\) and \(L_x=L_y=L\)
def solver(N, L, Nt, cfl=0.5, c=1, store_data=10,
u0=lambda x, y: np.exp(-40*((x-0.6)**2+(y-0.5)**2))):
xij, yij = mesh2D(N, N, L, L)
Unp1, Un, Unm1 = np.zeros((3, N+1, N+1))
Unm1[:] = u0(xij, yij)
dx = L / N
D = D2(N)/dx**2
dt = cfl*dx/c
Un[:] = Unm1[:] + 0.5*(c*dt)**2*(D @ Un + Un @ D.T)
plotdata = {0: Unm1.copy()}
if store_data == 1:
plotdata[1] = Un.copy()
for n in range(1, Nt):
Unp1[:] = 2*Un - Unm1 + (c*dt)**2*(D @ Un + Un @ D.T)
# Set boundary conditions
Unp1[0] = 0
Unp1[-1] = 0
Unp1[:, -1] = 0
Unp1[:, 0] = 0
# Swap solutions
Unm1[:] = Un
Un[:] = Unp1
if n % store_data == 0:
plotdata[n] = Unm1.copy() # Unm1 is now swapped to Un
return xij, yij, plotdata
Note that boundary conditions can be set after stepping the solution forward in time. Hence it does not really matter what is done to the solution at boundaries when Unp1 is computed, because the boundary will be modified later. This is very different from when we solve the Poisson equation and it is possible because of the recursive (explicit) solution algorithm. Note that the unknown Unp1 is not part of the right hand side in the update
A vectorized algorithm for this problem is possible, and it would read
However, this approach is not necessary, since we use a recursive solution algorithm that simply updates \(U^{n+1}\). Note that if we were to write this vectorized equation as Ax=b, then A would simply be equal to the identity matrix and b would be equal to the entire right hand side above. Since A is the identity matrix, the linear algebra problem is already taken care of. For the Poisson equation it is different because A equals \(D^{(2)}_x \otimes I_y + I_x \otimes D^{(2)}_y\). As stated initially, the hyperbolic nature of the wave equation makes it significantly easier to solve than the elliptic Poisson equation.
We are now ready to solve the wave equation for a chosen mesh size and we also choose to store the solution every fifth time step.
xij, yij, data = solver(40, 1, 501, cfl=0.5, store_data=5)
Plot the initial data
from matplotlib import cm
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
surf = ax.plot_surface(xij, yij, data[0], cmap=cm.coolwarm,
linewidth=0, antialiased=False)
Create an animation for the solution using matplotlib.
%%capture
# capture, otherwise there will be a plot in this cell
import matplotlib.animation as animation
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
frames = []
for n, val in data.items():
frame = ax.plot_wireframe(xij, yij, val, rstride=2, cstride=2);
#frame = ax.plot_surface(xij, yij, val, vmin=-0.5*data[0].max(),
# vmax=data[0].max(), cmap=cm.coolwarm,
# linewidth=0, antialiased=False)
frames.append([frame])
ani = animation.ArtistAnimation(fig, frames, interval=400, blit=True,
repeat_delay=1000)
ani.save('wavemovie2d.apng', writer='pillow', fps=5) # This animated png opens in a browser
from IPython.display import HTML
from IPython.display import display
display(HTML(ani.to_jshtml()))